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1 Introduction 



The process controlling the diferentiation of stem, or progenitor, cells into one specific functional direction is 
called lineage specification. An important characteristic of this process is the multi-lineage priming, which re- 
quires the simultaneous expression of lineage-specific genes. Prior to commitment to a certain lineage, it has 
O ' been observed that these genes exhibit intermediate values of their expression levels. Multi-lineage differen- 

r^ ' tiation has been reported for various progenitor cells |Hu et al., 1997[ |Akashi et al., 20031 |Kim et al., 2005| 

O^' [Miyamoto et al., 2002[ [Swiers et al., 2006[ [Loose fc Patient, 2006[ [Patient et al., 2007[ [Graf, 2002| , and it 

has been explained through the bifurcation of a metastable state [Roeder fc Glauche, 2006[[Huang et al., 2007[ 
[Chickarmane et al., 2009] . During the differentiation process the dynamics of the core regulatory network 
tH- . follows a bifurcation, where the metastable state, corresponding to the progenitor cell, is destabilized and 

'^ I the system is forced to choose between the possible developmental alternatives. While this approach gives 

fvj I a reasonable interpretation of the cell fate decision process, it fails to explain the multi-lineage priming 

^^ ' characteristic. Here, we describe a new multi-dimensional switch-like model that captures both the process 

^^ ' of cell fate decision and the phenomenon of multi-lineage priming. We show that in the symmetrical interac- 

tion case, the system exhibits a new type of degenerate bifurcation, characterized by a critical hyperplane, 
containing an infinite number of critical steady states. This critical hyperplane may be interpreted as the 
l^ . support for the multi-lineage priming states of the progenitor. Also, the cell fate decision (the multi-stability 

}_( , and switching behavior) can be explained by a symmetry breaking in the parameter space of this critical hy- 



perplane. These analytical results are confirmed by Monte-Carlo simulations of the corresponding chemical 
master equations. 

2 Stem Cell Differentiation 

The processes describing the interactions in systems like transcriptional regulatory networks are extremely 
complex. Genes can be turned on or off by the binding of proteins to regulatory sites on the genome 
[Ozbundak et al., 2002[ [Ptashne fc Gann, 2002| . The proteins are known as transcription factors, while the 
DNA-binding sites are known as promoters. Transcription factors can regulate the production of other tran- 
scription factors, or they can regulate their own production. The transcription process can be described by a 
sequence of reactions, in which RNA polymerase (R) binds to a gene's promoter leading to the transcription 
of a complete messenger RNA molecule. The genetic information transcribed into messenger RNA molecules 
is then translated into proteins by ribosomes. Thus, the general assumption is that the genes can be excited or 
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Figure 1 : The architecture of the self-excitation and inhibition mechanisms in the binary cell fate decision 
circuit. 

inhibited by the products of the other genes in the network, generating complex behavior like multi-stability 
and switching between different steady state attractors. Based on these general assumptions, it has been 
shown that a simple gene regulatory circuit (Fig. [1]) in which two transcription factors, X and Y , inhibit each 
other, and in the same time activate themselves, can be used as a model of binary cell fate decision in mul- 
tipotent stem or progenitor cells [Roeder fc Glauche, 2006[ [Huang et al., 20071 |Chickarmane et al., 2009] , 
This circuit can generate multistability and explains the symmetric precursor state, in which both factors 
are present in the cell at equal (low) amounts. This circuit typically produces three stable attractor states 
that correspond to observable cell states. The state 1, with the expression pattern X ^ Y, and the state 
2, with the opposite pattern Y :^ X represent the cell fates, while the state 3, with a balanced expression 
X ~ y represents the undecided multipotent state. This simple model provides a conceptual framework for 
understanding cell fate decisions, and it will be used as a starting point in the development of our model. 

3 Monte- Carlo Simulation Approach 

The Monte-Carlo simulation approach employed here is based on the well known Gillespie algorithm [Gillespie, 1977| , 
which is a variety of a dynamic Monte Carlo method. The traditional continuous and deterministic descrip- 
tion of biochemical rate equations, modeled as a set of coupled ordinary differential equations, relies on 
bulk reactions that require the interactions of millions of molecules. In contrast, the Gillespie stochastic 
algorithm simulates every reaction explicitly, and calculates the time evolution of the system by deter- 
mining the probabilities of each discrete chemical reaction and the resulting changes in the number of 
each molecular species presented in the system. This algorithm has rigorous theoretical foundations, and 
gives the exact solution for a system of elementary chemical reactions in the approximation of a well- 
mixed environment. When simulated, a Gillespie realization represents a random walk that exactly rep- 
resents the distribution of the chemical master equation. The algorithm is computationally expensive and 
several modifications have been proposed to speed up computation, including the next reaction method, 
tau-leaping, as well as hybrid techniques where abundant reactants are modeled with deterministic behav- 
ior [Gibson fc Bruck, 2000[ [Rathinam et al., 2003[ [Slepoy et al., 2008] . These adapted techniques provide 
a compromise between computational speed and the exactitude of the theory behind the algorithm as it 
connects to the chemical master equation. Here we use the standard stochastic simulation algorithm, known 
as the Gillespie's direct method. The rigorous derivation of the algorithm has been given elsewhere and it 
has been shown to remain "exact" for arbitrary low number of molecules [Gillespie, 1977 . 



Consider a system composed of N chemical species X^ (t^ = 1, ..., A^), interacting through M reactions 
i?p {fi — 1, ...,M) in the cell volume V. Every chemical reaction i?^ is characterized by its stochastic rate 
constant fc^, which depends on the physical properties of the molecules taking part in the reaction. The 
product k^dt is the probability that one elementary reaction R^ happens in the next infinitesimal time 



interval dt. The main steps of the Gillespie algorithm consist of: 

(a) calculating the waiting time t for the next reaction to occur; 

(b) determining which reaction fi in the system actually will occur. 

These quantities are computed by generating two random numbers according to the following probability 
density function: 

P{t, h) = Ufj, exp(-aoT), (1) 

where 

Qfj, = m^/cp, (2) 

and 

M 

ao^^afj.. (3) 

Here, m^ is the number of distinct reactant combinations available for the reaction i?^ at the given state of 
the system. The coefficient a^ is called the propensity of reaction R^. Thus, P{t,ijl) is the probability that 
the next reaction will occur in the infinitesimal time interval dt and that it will be the R^ reaction. After 
determination of r and /i, the numbers of molecules in the system are adjusted according to the reaction 
i?^. Also, the time t is advanced to i + r. The larger the propensity is, the greater is the chance that a 
given reaction will happen in the next step of the simulation. It is worth noting that there is no constant 
length for a time-step in the simulation. The length of each time-step is determined independently in every 
iteration, and takes different values depending on the state of the system. 

The implementation of the Gillespie algorithm is straightforward, and one can find excellent descriptions 
of it in the literature [Adalsteinsson et al., 2004} [Kierzek, 2002| . Below we give the pseudo-code of the 
algorithm: 

T^Gillespie's direct method 

1. Set initial numbers of molecules, set time t ^ 0; 

2. Galculate the propensities, a^, for all /i = 1, ...,M; 

3. Choose /i with the probability: 

Pr{reaction = /i) = — jj^ ; (4) 

4. Choose T with the probability: 



Pr(iime = t) = 2^ ^, 



M 

"v- 



exp 




(5) 



5. Change the number of molecules to reflect execution of reaction /i; 

6. Set i <— i -I- T, and go to step 2. 

4 2-Diniensional Model 

We consider the two gene circuit shown in Figure 1. We will focus on the elementary processes that must 
occur, such as the promoter binding of the transcription factors X and Y to the promoters, A and _B, re- 
spectively, and the activation and degradation of transcription factors. Also, we propose a general approach 
to integrate the two inputs to each gene, which does not depend on the assumption of cooperativity or other 



explicit modeling. In order to provide a quantitative model of this genetic circuit, we employ a formalism 
originally developed for the mean-field description of the stochastic interactions in transcriptional regula- 
tory networks [Andrecut fc Kauffman, 2006[ [Andrecut et al., 2008| . The promoter binding and unbinding, 
subsequent self-activation, inhibition, dissociation and the degradation reactions for X, and respectively Y, 
are: 



k 



AR 



(6) 
A + R — > A + R + X 



'^AX 



(7) 

A + X ^f AX 



^AX (8) 

AX ^ A + X 



fet 



(9) 



AX + R — > AX + R + X 

(10) 



k+ 



A ^ Y — > AY 



AY — > A + Y 



k- 



X 



ksR 

B + R^B + R + Y 



k+ 
B + Y — > BY 



kgY 

BY — > B + Y 



BY + R — > BY + R + Y 

k+ 

^BX 

B + X — > BX 
BX — > B + X 



(11) 
(12) 
(13) 
(14) 
(15) 
(16) 
(17) 
(18) 
(19) 



y ^ 

Here, kj^x' ^AX' ^bY' ^by^ describe the binding and release rates between the transcription factor and 
the promoter element, fcj^y, fc^y, ^bX' ^bx correspond to the cross inhibition rates, while fcj, k^, ky, ky 
reflect the activation and the degradation rates of the transcription factors. We assume that the role of the 
first reaction for each transcription factor. Equation [S] and respectively Equation [T^l is just to provide a 
small "basic level of expression" (with the rates kAR, and respectively ksR), in order to avoid their complete 



extinction. It's effect is equivalent to a positive noise term tjx.y in the differential equation describing the 
dynamics of the transcription factor. Therefore, in the following analysis we will neglect the contribution of 
this reaction, since it doesn't really have an influence on the "logic functionality" of the circuit. 

The dynamical behavior (rate of change of active levels of the proteins) of the isolated transcription 
factors is therefore described by the stochastic differential equations: 



d 
It 



-[X] = k+[AX]-k-^[X]+i^x. (20) 



j^[Y]^k+[BY]-ky[Y]+^y, (21) 

where [.] denotes concentration. Assuming that the reversible binding-unbinding processes are in equilibrium, 
we have: 

kAx[A][X] = [AX], (22) 

kBY[B][Y] = [BY], (23) 

kAY[A][Y]^[AY], (24) 

kBx[B][X] = [BX], (25) 

where kAx — k^x/^AX^ ^by — kgy/^BY^ ^ay = ^ay/^\y^ ^bx — kgx/^BX- Also, since the promoters 
can be in three different states we have: 

[AX] + [AY] + [A] = [Ao], (26) 

[BY] + [BX] + [B] = [Bo], (27) 

where [Aq] and [Bq] are the total concentrations of the two promoters. From the above equations, and 
neglecting the noise terms, we obtain the following system of deterministic differential equations: 

l:c^ax( ^ --iV (28) 

at \aix + a^y + 1 / 

at \bix + h2y+l J 

where we assumed that: x = [X], y = \Y], a = k^, j3 = ky, ai = kAx, &i = ksx, 12 — kAY, &2 — ksY, 

as = [AQ]kAxkx/kx, &3 == [BolfcsYfcy/fcy 

We are interested in the symmetrical case, where a = /3, ai = &2, 02 = &i, 03 = &3, such that the system 

becomes: 

—X — ax ( 1 I , (30) 

at \aix + a2y + 1 / 

±y = ay( ^ iV (31) 

at \a2X + aiy + 1 / 

The steady states of the above differential system of equations are given by the solutions of the non-linear 
system: 

^x = 0^ Fix, y, a, {a}) = ax ( -^^—7 - l) = 0, (32) 

at \aix + a2y + 1 / 

^y = ^ G{x, y, a, {a}) - ay ( -^^— - - l) = C (33) 

at \a2X + aiy + L J 



In this case, one can easily verify that the system has four steady states: 



(a:o,2/o) = (0,0) 
'as - 1 



[xi,yi) 



ix2,y2) = 0, 






(2^3,2/3) 



(03 - 1,03 - 1), 



(34) 
(35) 
(36) 

(37) 



ai + 02 

corresponding to the extinction, exclusive and coexistence equilibria. These fixed points are positively defined 
if as > 1. 

In order to evaluate the local stability we calculate the eigenvalues, A and /i, of the Jacobian matrix at 
these steady states: 



J{x,y,a,{a}) = 



dF OF 

dx dy 

dG dG 

dx dy 



a 



{aix+a-2y+iy 



- 1 



{aix+a2y+l)'^ 
a3{a2X+l) _ -, 



{a2X+aiy+lY (02^+01^+1) 

The eigenvalues of the Jacobian for the extinction state (xo, j/o) are: 

A = a{as - 1) > 0, 
^ ~ a{as — 1) > 0. 



(38) 



(39) 
(40) 



Thus, this steady state is always unstable, since 03 > 1. The eigenvalues for the exclusive steady states, 
(xi,2/i) and (2:2,2/2), are: 



03-1 



and respectively: 



H = —a 



A = —a 



03 




(03 - l)(a2 - 


' Oi) 


Oi + 02(03 - 


-1) ' 


(03 - l)(a2 - 


-oi) 


ai + 02(03 - 


-1) ' 


03 — 1 

— — rv _ 





03 



(41) 
(42) 

(43) 

(44) 



Therefore, the exclusive equilibria are stable if 02 > oi, and unstable if 02 < Oi. In contrast, the eigenvalues 
for the coexistence equilibrium (0:3, 2/3) are: 



A = —a 



fi = —a- 



(03 - 1)(Qi +Q2) 
03(01+02) 

(03 - l)(ai - 02) 



(45) 



(46) 



03(01 +02) 

Since A < 0, this steady state is stable ii fi < 0, and it looses stability if /i > 0. One can easily see that 
the stability condition, /i < 0, is equivalent to 02 < oi. Thus, a change in the ratio p ~ 01/02, triggers a 
bifurcation from one stable steady state (xs, ys), when p < 1, to two stable steady states (xi, yi) and (2:2, 2/2)) 
when p > 1. 

In Fig. [2]and Fig. [3]we give the results of the Monte-Carlo simulations. The initial concentrations and the 



main reaction constants are set as: B = 100, Aq = Bq = 1, Xq — Yq = Q, k^ = ky = 0.01, kA = ks = 0.01, 
^x = ky = 0.01, k'^^ = k^Y = Ij k^Y ~ ^BX ~ 1- ^^S- m gives the trajectories x{t) — X{t)/B and 
y{t) = Y{t)/ B for /? < 1, when there is one noisy attractor, corresponding to the coexistence equiUbrium, 
(2^3, 2/3) (Fig. [UJa)), and for p > 1, when there are two noisy attractors, corresponding to the exclusive 
equilibria, (a;i,2/i) and (2:2,2/2) (Fig. HKb)). Also, in Fig. [3] we have represented graphically the probability 
density distribution, P{x,y), of the transcription factors (obtained by averaging over Af = 10'' trajectories 
with T = 10^ reactions events). One can see that for p < 1, the system has only one noisy attractor, 
corresponding to the stable fixed point (2:3, j/s) (Fig. [Sj^a), ai — 1, a2 — 2), while for p > 1, the system 
exhibits two noisy attractors corresponding to the stable fixed points (xi,yi), and respectively {x2,y2) (Fig. 
^h), fli = 2, 02 = 1). We should note that the absolute values of the rate constants do not play a critical 
role in the simulation, as long as their ratios satisfy the bifurcation constraints. 

An important case of the above analysis corresponds to the critical bifurcation parameter p = 1. In this 
case the system has the form: 

—x = ax^{x,y,{a}), (47) 

—y^ay^{x,y,{a}), (48) 



where 



Hx,y,{a})^ l\ 1. (49) 

ai[x + y) + 1 



One can easily verify that in this case, the exclusive and coexistence equilibria disappear, and the system 
has an infinite number of stedy states fi = {{x, y) G R^|$(2;, y, {a}) = 0}, which are practically equivalent 
to the positive segment of the linear equation: x + y — [a^ — l)/a\. These steady states have the following 
eigenvalues: 

A = 0, (50) 

p = -(03 - l)/a3 < 0. (51) 

Therefore, the steady states Q, are stable, and the system undergoes a degenerate bifurcation (see Appendix). 
This situation is presented in Fig. [SJc) and Fig. [HJc), for ai = 02 == 1. One can see that the stochastic system 
is "undecided", exploring every point of the critical line with non-zero probability. The line is attracting, 
except along itself, that is, there is no "longitudinal" force on this line. Therefore every state on it is 
indifferently stable. Thus, the critical line becomes an ergodic attractor. 

5 AT-Dimensional Model 

We consider a A^— gene circuit, where we denote by A„ the transcription factors and by yl„ the promoters, 
n = 1, ..., A^. For each gene we assume the following set of equations: 

'=" (52) 



k+ 



A.„ + A„ 






(53) 
(54) 





V. 0.5 






(C) 



Figure 2: Monte-Carlo simulation of the 2-dimensional circuit: (a) p < 1, one attractor {k^x ~ ^by ~ ^-^ 
and k~^Y — ^bx ~ 1)' (^) P > 1; two attractors {k^x ^ ^by ^ ^ ^^^'^ ^ay ^ ^bx ^ 0-5); (c) p = 1, the 
critical case of the degenerate bifurcation {k^x — ^^by ~ ^ay ~ ^bx ~ ^)- 






Figure 3: The probability density distribution, P{x,y), of the stochastic trajectory of the 2-dimensional 
circuit: (a) p < 1, one attractor; (b) p > 1, two attractors; (c) p = 1, the critical case of the degenerate 
bifurcation. 



'=" (55) 

AnXn + R > AnXn + R + Xn 

''" . (56) 

Thus, the promoter A„ can bind to any of the N transcription factors. Thus, the dynamical behavior of the 
transcription factors is described by the foUowing system of stochastic differential equations: 

^^[X„]=k+[A„Xn]-k-[Xr,]+vx„, n = l,...,iV, (57) 

where ?7x„ is the noise term corresponding to the first reaction (51). Assuming that the reversible binding- 
unbinding processes are in equilibrium, we have: 

knm[An][X„,] ^ [A„X„], n,m^l,...,N, (58) 

where k„m = ^nml^nm- Also, since the promoters can be in A^ + 1 different states we have: 

N 

Y, [AnX^] + [An] = K], n = 1, ..., M. (59) 

m— 1 

where [A°] are the total concentrations of the promoters. From the above equations, and neglecting the 
noise terms, we obtain the following system of deterministic differential equations: 

d I Pnk., 



dt 



Xn^anXn\—-T:j 1 , n = l,...,N. (60) 



A^?n— 1 '^nni'^'m ' ^ 



where we assumed that: x = [Xn], ctn — fc„ , /3„ — [^°]fc^/fc„ • 
Let us consider now the symmetric case: 

d ( Bk \ 

— a;„ = F„(xi,...,xjv,a,/3, «;,7) = aa;„ ■ ^^ — - - 1 , n = l,...,iV, (61) 

dt V«2;„+7^„_,„x„ + l j 

where k ~ knn and 7 = knm for m ^ n = 1,...,N. The steady states corresponds to the solutions of the 
nonlinear system: 

— a;„ = 0^ F„(xi,...,X7v,a,/3, K,7) = 0, n=l,...,iV. (62) 

There are N + 2 steady states: 

(:c(°\...,xi^))=(0,...,0), (63) 

(x«, ..., x«, ..., x«) = (^0, ..., ^^, ...0) , ^^l, ..., N, (64) 

Again, these states correspond to extinction, A^— exclusive and coexisting equilibria, and they are positively 
defined if /3k > 1. The stability of these states can be analyzed using the eigenvalues of the Jacobian matrix: 



J(xi,...,xjv,a,/3, K,7) 



d 

—Fn{xi,...,XN,a,(3,K,j) 



(66) 

■.A=1....,N 
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where 

7 ^ Y? ~ Q; ij n — I 

('^^,.+7E„^„^r„+ij , (67) 

_ ^Pkjx^ ^j ^^■ 

The eigenvalues of the extinction state are: 

A„ = a(/3K-1), n^l,...,N, (68) 

which means that this steady state is always unstable, since /3k > 1. For the exclusive equilibria the 
eigenvalues are: 

K = \ .(,.-^i1(,-.) / , n=l,...,iV. (69) 

Therefore, these states become stable if 7 > k, and unstable if 7 < k. In the case of coexisting equilibrium 
the Jacobian is given by: 




a(l-ffK) 

l'^+(A^-i; 
Q7(/3k — 

/jfc[K+(Ar-i)7] 



and it has the following eigenvalues: 







n=l,...,N. (71) 



a(ffK-l)(K-7) . ^ „ 3^ 7 



Thus, this state becomes stable if 7 < k, and unstable if 7 > k. 

In the critical case, k — j, the steady state equations are degenerated and we have again an infinite 
number of steady states, all of them satisfying the critical hyperplane equation: 



In this critical case the Jacobian takes the simplified form: 

dxi " li 

and it has the following eigenvalues: 



o -^n — /s*^"' v'"^/ 



^^ ., ,. , n = l,...,iV. (74) 

ij n ^ I 

Thus, one eigenvalue is always negative, since /3k > 1, and the other iV — 1 eigenvalues are zero. Therefore, 
the hyperplane containing the infinite number of steady states is attractive and marginally stable. 

In Fig. [5] we give the simulation results for a circuit consisting of three genes, A^ = 3. The initial 
concentrations are set as i? = 150, Ai = A2 = A3 — 1, Xi — X2 ~ X^ — 0. The rate constants are the 
same as for the 2-dimensional circuit. Fig. |4] gives the trajectories Xn{t) for 7 < k, when there is one noisy 
attractor, corresponding to the coexistence equilibrium, {x\ ,X2 ,2:3 ) (Fig. Hl^a)), and for 7 > k, when 
there are three noisy attractors, corresponding to the exclusive equilibria, {x^ , 0:2 i Xg ) = I ^^ , 0, J , 

(xf\x^^\4^^) = (O'^ir^'O) and (xf\x^^\x^^^) = ^0,0, ^^] (Fig. g^b)), and in the degenerate case 
when the plane xi + X2 + x^ ~ ' is the ergodic attractor (Fig. IHJc)). 
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Figure 4: Monte-Carlo simulation of the 3-dimensional circuit: (a) 7 < k, one attractor; (b) 7 > k. three 
attractors; (c) 7 = k, the critical case of the degenerate bifurcation. 

6 Conclusion 

We have presented a new inulti-diinensional switch-like model that captures both the process of cell fate deci- 
sion and the phenomenon of multi-lineage priming. The previous attempts to model the coexistence of three 
discrete stable states are based on a complicated molecular interaction mechanisms, requiring cooperativity 
or additional transcription factors |Roeder fc Glauche, 2006| [Huang et al., 20071 |Chickarmane et al., 2009] . 
Here, we have shown that very elementary cross-inhibition between two genes and independent autoactiva- 
tion can give rise to multistability without cooperativity. It is important to note the obvious fact that the real 
molecular mechanisms that govern the dynamics of this gene regulatory circuit are by orders of magnitudes 
more complex, involving perhaps thousands of steps not accounted for in the presented model. However, 
this simplified description is still able to capture to qualitative cell fate decision behavior, specifically, the 
existence of an indeterminate multi-potent progenitor state with equal levels of transcription factors, and 
the generation of stable attractor states with asymmetric expression patterns. Also, we have shown that in 
the symmetrical interaction case, the system exhibits a new type of degenerate bifurcation, characterized by 
a critical hyperplane containing an infinite number of critical steady states. This degeneration of the central 
attractor state captures the intrinsic heterogeneity of the undecided multipotent state allowing individual 
cells to a range of states, and may be interpreted as the support for the multi-lineage priming states of the 
progenitor. Also, the cell fate decision (the multi-stability and switching behavior) can be explained by a 
symmetry breaking in the parameter space of this critical hyperplane. It is important to note here that the 
critical hyperplane is also ergodic. Thus, in the critical regime, any stochastic trajectory of the system will 
be attracted to the critical hyperplane. Also, in this case, the dynamics will become confined to this region, 
such that the system will visit all the points of the critical hyperplane with non-zero probability (the priming 
phenomenon). However, any perturbation of this critical hyperplane will force the system to collapse in one 
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of its non-trivial stable steady states (the cell fate decision process) . 

7 Appendix: Degenerate Steady State Bifurcation 

We consider a 2-dimensional system of stochastic differential equations (SDE) , of the following generic form: 

\ y = /3(y-y)*(a;,y, {7})g(a;,y, {6}) + ?7y 

where $,/, 5 : M^ — > R , x, j/,x,y, a,/3, {7}, {a}, {&} G R, and rj^, rjy correspond to additive noise terms. 
Also, we denote hy il — {{x,y) G R^|<i>(x,2/, {7}) = 0} the set of solutions of the equation $ = 0, and we 
assume that: f{x, y, {a}) ^ and g{x, y, {b}) ^ 0, for any (a;, y) G R^. 

Theorem 1: The SDE system {S) exhibits a degenerate bifurcation, (x, y) — ;> i7 C il, from one steady 
state (x, y) to a subset of steady states fi C Jl, if: 

max{-a$(x, y, {7})/(S, y, {a}), -/?$(x, y, {l})g(x, y, {6}} > 0, (76) 

and 



for any (x, y) G fi and 



V,„$(x, y, {7}) = («, V$(x, y, {7})) < 0, (77) 



V = [a(x - x)/(x, y, {a}), /3(y - y)g(x, y, {6})]'^. (78) 



Proof: A first steady state of the system (S) is (x, y). Also the system has an infinite number of steady 
states 51, corresponding to the solutions of the equation $(x,y, {7}) — 0. The stability of the steady states 
can be analyzed using the eigenvalues, Aq and Ai, of the Jacobian matrix J, which are given by the solutions 
of the equation | J— A/| = 0, where / is the identity matrix. In general, the eigenvalues are complex numbers, 
and the distance between the solution of the system and thesteady state changes at an exponential rate, given 
by the real part of the eigenvalue. For simplicity the following discussion is restricted to real eigenvalues, 
though steady states with complex eigenvalues have similar properties based on the value of the real part 
of the eigenvalue. A negative eigenvalue implies that the solution approaches the steady state along the 
corresponding eigenvector, while a positive eigenvalue implies that the solution moves away from the steady 
state along the eigenvector. In a 2-dimensional system there are three possible cases. A stable steady state 
has two negative eigenvalues, and hence attracts all the solutions in a surrounding region. An unstable 
steady state has two positive eigenvalues and all the solutions in its neighborhood move away from it. A 
saddle point has one negative and one positive eigenvalue. Now let us analyze the stability of the steady 
states of the system. The eigenvalues of the Jacobian for (x, y) are: 

I Ao = -a/(x,y, {a})$(x,y, {7}) 
\ Ai = -/35(x,y,{&})$(x,y,{7}) 

Thus, this steady state is stable if Ao, Ai < 0, and it looses stability if Aq > or Ai > 0, which is equivalent to 
the condition imposed by the EquationlTSl The other steady states (x, y) G 51, have the following eigenvalues: 



Ao=0 

•a* I _ 

Qx ~^ f-'\y i/yy Qy 



Ai^a(x-x)/||-H/3(y-y)g-f ' ^^"^ 



and they are degenerated, since they have at least one zero eigenvalue. Also, these degenerate steady states 
become stable for Ai < 0. Since Ai — {v, V<i>) — Vi,$, this stability condition is equivalent to the condition 
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imposed by the Equation [771 which requires that the derivative of $ in direction v must be negative. The 
directional derivative of a manifold $ along a vector w at a given point {x,y), intuitively represents the 
instantaneous rate of change of the manifold, moving through {x,y), in the direction of v. One can easily 
verify that in this case, v is the eigenvector of the Jacobian corresponding to the eigenvalue Ai, that is we 
have: Jv = Xiv. Thus, any change in the parameters {7}, such that the stable steady state (x,y) becomes 
unstable, and the steady states fl = {{x,y) £ ft\\7y^{x,y,{j}) < 0} become stable, results in a degenerate 
bifurcation of the dynamics of the stochastic system (S) . 

A similar property can be formulated for stochastic discrete maps (SDM) of the following generic form: 

a<^{xt,yt, {l})f{xt,yt, {a})] + x + t]^ 
I3<^{xt,yt, {'i})9{xt, yt, {b})] + y + Vy 

where $, /, 5 ; R^ — )■ R , xt,yt,x,y,a, l3,{j}, {a},{b} G M, and r]^, rjy correspond to additive noise terms. 
We denote by J7 = {{x, y) e M^|$(a;, y, {7}) — 0} the set of solutions of the equation $ = 0. Also, we assume 
that: /(x, y, {a}) / and g{x, y, {b}) ^ 0, for any {x, y) G M^. 

Theorem 2: The SDM system (Af) exhibits a degenerate steady state bifurcation, {x,y) — > fi C il, 
from one steady state (a;, y) to a subset of steady states il C fi, if: 

max {|1 - a$(S, y, {^])f{x, y, {a])\ , |1 - /3$(x, y, {i})g{x, y, {h])\] > 1, (82) 

and 

|l-V„$(x,2/,{7})|<l, (83) 

for any (x, y) G 51 and 

V = [a{x - x)fix, y, {a}), (3{y - y)g{x, y, {b})f. (84) 

Proof: The steady states of the map (P) are (S, y), and the set 51, corresponding to the solutions of the 

equation $(x, y, {7}) = 0. As before, the eigenvalues, Aq and Ai, are given by the equation \J — A/| = 0, 

where J is the Jacobian of the discrete map (P). In the case of discrete maps, a stable steady state is 

characterized by |Ao| < 1 and |Ai| < 1, while an unstable steady state is characterized by |Ao| > 1 or 

All > 1. The eigenvalues of the Jacobian for (x,y) are: 

a$(x,y,{7})/(x,y,{a}) 
-/3$(x,y,{7})5(x,y,{6}) 

Thus, this steady state is stable if |Ao| < 1 and |Ai| < 1, and it becomes unstable if |Ao| > 1 or |Ai| > 1, which 
is equivalent to the condition imposed by the Equation 1821 The other steady states 11, have the following 
eigenvalues: 

r Ao-i 

\ Ai = l-a(x-i)/|f-/3(y-y)yf 

and they are degenerated, since they have at least one eigenvalue equal to one. These degenerate steady 
states become stable for |Ai| < 1, which is equivalent to the condition imposed by the Equation l83l Also, 
one can verify that v is the eigenvector of the Jacobian, corresponding to the eigenvalue Ai. Thus, any 
change in the parameters {7}, such that the stable steady state (x,y) becomes unstable, and the steady 
states 51 = {(x, y) G 51| |1 — V«$(a;,y, {7})! < 1} become stable, results in a degenerate bifurcation of the 
dynamics of the stochastic discrete map (M) . 




a$ n, -X r)* ' (86) 
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